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We evaluate perturbatively the density matrix in the low-temperature limit and thus the ground- 
state wave function of the anharmonic oscillator up to second order in the coupling constant. We 
then employ Kleinert's variational perturbation theory to determine the ground-state wave function 
for all coupling strengths. 
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I. INTRODUCTION 

Variational perturbation theory as developed by Kleinert JjJ provides a systematic algorithm to evaluate pertur- 
bation series at all coupling strengths including the strong-coupling limit g — ► 00. It was thoroughly investigated for 
the ground-state energy of the anharmonic oscillator up to the 150th order Bpl, and its convergence was found to 
be exponentially fast and uniform. A similar systematic study has notyet been performed for the ground-state wave 
function. A first-order variational approach was set up by Kunihiro ]4j. However, his method did not satisfactorily 
deal with certain problems in the variational procedure which will be discussed in detail below. 

In this work we improve Kunihiro's first-order calculation and extend the treatment to the second order in the 
coupling strength. The ground-state wave function of the anharmonic oscillator is calculated from the low-temperature 
limit of the diagonal elements of the density matrix. A variational evaluation of the density matrix for the double-well 
potential has already been perfomed for finite temperatures in Ref. j5j. 

II. PERTURBATION THEORY 

Consider a quantum-mechanical point particle of mass M moving in the one-dimensional anharmonic oscillator 
potential 



(2.1) 



where uj denotes the frequency and g the coupling constant. We determine its ground-state wave function ^(x) by 
evaluating the low-temperature limit of the diagonal elements of the density matrix 



ty(x) = lira \J p(x, x) , 

j3—>oo 



which is defined by 



p(x b ,x a ) 



(x b hj3\x a 0) 



x(hf3)=x b ( ^ ph0 

T>x exp < — — / dr 

x(0)=x a I "JO 



^(r) + ^¥(r)+^(r) 



Here (xf, h[3\x a 0) denotes the imaginary-time evolution amplitude with the path integral representation |jj 

(x b h(3\x a 0) 

and Z denotes the partition function 

/ + OC 
dx (xhf3\x0) 
-00 



By expanding Eq. (2.4) in powers of the coupling constant g we obtain the perturbation series 

i-h/3 „2 phl3 r hf3 



(x b Hf3\x a 0) = (x b h(i\x a 0) u 



1_£ / dnixHr^ + £2 / dT i / ^(.x 4 (Ti)x 4 (r 2 )), 
h Jo 2h Jo Jo 



(2.2) 



(2.3) 



(2.4) 



(2.5) 
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where we have introduced the harmonic imaginary-time evolution amplitude 



(xb fif5\x a 0) 



x(h/3) = x b ( ^ rhp 

Vx cxp < — — / dr 

x(0)=x a I "JO 



M . M 2 2 . . 
( T ) + y^ 2 (r) 



and the harmonic expectation value for an arbitrary functional F[x] of the path x(t): 



(x b Ti(3\x a Qi) u J x (o) = 



x(hf3)—Xb 



VxF[x] exp < — \ 



Tift 



dr 



M . 2 . . M 2 2 , . 
yX 2 (r) + - W V(r) 



(2.7) 



(2.8) 



The latter is evaluated with the help of the generating functional for the harmonic oscillator, whose path integral 
representation reads 



(x b tiP\x a 0) u \j] = 



x(h/3)=Xb 



x(0)=x a 



Tift 



dr 



M -2/ \ M o 2l s ■ / \ / \ 



(2.9) 



leading to 

{x b h(i\x a $)u\j] = (xb h/3\x a 0) UJ ex.p 
with 

(x b h(3\x a 0) 



I rhft 1 ,hft M3 

ft Jo dTlXcl<yTl ^^ + 2tfJ dTl J dT2 G ( Tl ' T2 )-7'( Tl )i( r 2/ 



Mlj 



exp 



27r7i sinh U(3u) \ 2hsmhfi(3ui 
In Eq. ( p.lC ) we have introduced the classical path 

x a sinh(?i/3 — r)a; + x;, sinh ut 



[(x 2 + x 2 ) cosh K/3lu — 2x a Xb] 



sinh h/3u> 



and the Green function 
G(n,r 2 ) '' ' 



2M uj sinh Ti/Jw 



[0(ti — t 2 ) sinh(?i/3 — ti)u> sinhtJT2 + 9(t 2 — t\) sinh(?i/3 — t-i)uj sinhcjri] 



(2.10) 



(2-11) 



(2.12) 



(2.13) 



We fo llow Ref. || and evaluate harmonic expectation values of polynomials in x arising from the generating functional 
( [2.10 ) according to Wick's theorem. Let us illustrate the procedure to reduce the power of the polynomial by the 
example of the harmonic expectation value 



(x n (T 1 )x m (T 2 )) u 



(2.14) 



(i) Contracting x(n) with x n 1 (ri) and x m (T2) leads to a Green function G{t\,t\) and G(ti,T2) with multiplicity 
n — 1 and m, respectively. The rest of the polynomial remains within the harmonic expectation value, leading 
to {x n - 2 (n)x m (T 2 ))u J and (x n - 1 (Ti)* m ~ 1 fa)}<-- 

(ii) If n > 1, extract one x(t\) from the expectation value giving x c \(t\) multiplied by (x n ~ 1 (Ti)x m (T 2 )} UJ . 
(hi) Add the terms (i) and (ii). 

(iv) Repeat the previous steps until only products of expectation values (x(ti)) w = x c i{t±) remain. 
With the help of this procedure, the first-order harmonic expectation value (a; 4 (r 1 )) w is reduced to 

(x 4 (n)) w = zci(n) (^(n))^ + 3G(t 1 ,t 1 )(x 2 (t 1 )) uj . 

Furthermore, we find 

(a; 3 (ri)) w = x cX {t 1 )(x 2 {t 1 )) uj +2G(n,n)x c i(n), 



(2.15) 
(2.16) 
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and 



(x 2 (n)) u = x^(n) + G(n,n) . 



Combining Eqs. (2.15)— (2.17) we obtain in first order 

{x\ti)) u = 4(n) +6 Ic 2 1 (r 1 )G(T 1 ,T 1 ) + 3G 2 (r ll r 1 ). 
The second-order harmonic expectation value requires considerably more effort and finally leads to 



(2.17) 



(2.18) 



(x A (n)x 4 (r 2 )) u = a&nja&fo) + 16 x 3 c} (n) G(t u t 2 ) x 3 c} (t 2 ) + 12x 2 cl (n) G(r u n) 4(r 2 ) 

+72 x 2 cl (n ) G 2 (n , r 2 ) a^, (r 2 ) + 36 x\ (n ) G(n, n) G(r 2 , r 2 ) x\ (r 2 ) + 9 G 2 (n , n ) G 2 (r 2 , r 2 ) 
+96 x\ (n ) G(n , r 2 ) G(r 2 , r 2 ) z cl (t 2 ) + 6 G 2 (r x , n ) 4 (r 2 ) + 96 x cl (n ) G 3 ( n , r 2 ) x cl (r 2 ) 
+144z c i(Ti) G( n , n) G(n,r 2 ) G(r 2 , r 2 ) x cl (r 2 ) + 36G 2 (t!, n ) a&fa) G(r 2 , r 2 ) 
+144a: c 2 1 (T 1 )G 2 (T 1 ,T 2 )G(T 2 ,T 2 ) + 72G(T 1 ,T 1 )G 2 (T 1 ,r 2 )G(T 2 ,r 2 ) + 24G 4 (T 1 ,T 2 ). (2.19) 

The contractions can be illustrated by Feynman diagrams with the following rules. A vertex represents the integration 



/ dr, 

Jo 



(2.20) 



a line denotes the Green function 



G(ti,t 2 ) , 



(2.21) 



and a cross pictures a classical path 



1 = x c i(n) 



(2.22) 



Inserting the harmonic expectation values (218) and (2.19) into the perturbation expansion (2.£) leads in first order 
to the diagrams 



tip 



dTl(x 4 (Tl)) u 



+ 6 x- 



(2.23) 



whereas the second-order terms are 



dn I dT 2 {x i (Ti)x i (T 2 )) L j = 
Jo 



+36 x x x ^ - ^ x +96 x- 



+36 



144 





+ 96 x Q x +144 x- 




(2.24) 



We observe that contributions from both connected and disconnected Feynman diagrams appear. The disconnected 
diagrams vanish once we rewrite the imaginary-time evolution amplitude in the form 



{xb h(3\x a 0) = {xb hp\x a 0)uj exp [W(ai6, Tif3; x a ,0)} 



(2.25) 



where t he e xponent W(xb, hp; x a , 0) contains only the connected Feynman diagrams. We obtain from ( |2.6| ) and 
(|2.23Dd2T25|) the expansion 
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+48 x- 



-*+48 x- 



-x + 72 x- 



27^ 



36 




-x +72 




36 
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(2.26) 



where disconnected diagrams are indeed no longer present. As mentioned above, we restrict ourselves to the low- 
temperature limit of the diagonal elements of t he de nsity matrix which determin e the ground-state wave functi on. In 
order to evaluate the various contributions in (2.26), we need the classical path ( [2.12 ) and the Green function ( 2.13| ) 
in the low-temperature limit: 



lim x c i(r) = x (e^ T + e -*>W-r)) , 



lim G(n,T 2 ) = — — 
p—>ac 2Mui 



In - t 2 



) e -"(n-r 2 ) + 0( T2 _ n ) e - W (T 2 -n) _ e -^(r 1+ r 2 ) 



-2n/3w+w(n+T 2 ) 



(2.27) 

(2.28) 



Computing with thes e exp ressions the Feynman diag rams in (2.26), the low-temperature limit of the imaginary-time 
evolution amplitude ( 2.25| ) reads together with (2.11) 



lim (xhf3\xQ) = lim \i — — exp 

/3^oo 0^oo V rnr 



Muj 



hj3oj Muj 



9 ( 9K 2 



3h 



h \8M 2 uj 3 AM 2 lo 2 2Muj 2 ' 



2uj 



+ 



2057T 



2lh 6 



g*_ ( 21h 5 f3 
2h 2 V4Af 4 w 5 16Af 4 w 6 ' 2M 3 lo 5 ' 



nr 

4A/ 2 a 



r x 4 + 



3Afw 3 ' 



(2.29) 



According to (2.5), the partition function Z follows from (2.2E) by performing an integration with respect to x. This 
results in 



lim Z = lim exp 

13 — >oo — >oo 



Ti(3uj Zgti z (3 2lg 2 h 3 f3 



4M 2 lu 2 



8M 4 u 5 



(2.30) 



Inserting ( j2.29| ) and ( |2.30D into (^J) we observe a cancellation of all terms which would diverge in the low-temperature 
limit [3 — ► oo. Thus the diagonal elements of the density matrix read in this limit 



lim p(x, x) = J exp 
/3->oo V nir 



Muj 2 g ( 9h z 3h 



g 

2?i 2 



205?i 4 2ir 
16M 4 w 6 + 2A/ 3 c i 



Ilk 2 



1 

' - T — IE 

2w 
ft 



(2.31) 



By taking the square root and expanding the expo nent ial term up to second order in the coupling strength g, we 
derive the second-order ground-state wave function (fe^ 



T . , /Muj\ 1/4 ( Muj 

*^ = ("^rj ex p'-— 



2/T 



2h 

1559ft 4 141ft 3 



97T 



3ft 



16Af 2 w 3 4Afw 2 



1 

4o/ 



53ft 2 



256Af 4 w 6 32M 3 o 



32M 2 w 4 ' 



13ft 

24Afw 3 " 



16w 2 



(2.32) 



This result corresponds to the solution of the Bender- Wu recursion relation Q for the ground-state wave function 
which is normalized such that 



dx^ 2 (x) = 1 

-oo 

holds up to second order in the coupling strength g. 



(2.33) 
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III. VARIATIONAL PERTURBATION THEORY 



Variational perturbation theory enables us to evaluate the ground-state wave function for all values of the coupling 
constant g and even in the strong-coupling limit g — > oo. To this end we simply add and subtract a harmonic oscillator 
of trial frequency Q to the anharmonic oscillator potential (|2.lD : 



V(x) = 



M 

Y 



n 2 x 2 



m to 2 - n 2 

'~2 9 



gx 



(3.1) 



We now treat the second term as if it was of the order of the coupling constant g. T he result is obtained most simply 
by substituting the frequency uj in the original anharmonic oscillator potential (2J.) according to Kleinert's trick [Q 



gr, 



where we have 



n 2 



gW 



(3.2) 



(3.3) 



Writing the ground-state wave function (2.32) in the form ^(x) = exp[VF(x)] with the cumulant expansion 

= 2 



W{x) 



exp 



1 (Mw\ Muj 2 g ( 9h 
4 ° g \J!^~ J ~ ~2h X + h \16M 2 uj 3 4Muj 2 
2lh 4 , 



3h 2 1 t 

x ~l^ x 



2h 2 



2057T 



32M 4 w 6 4M 3 u 5 ' 



8MV 



6Moj 3 



(3.4) 



we a pply the trick (3.2) to W(x), and reexpand in powers of g at fixed r. Afterwards r is substituted according to 
( |3.3|) . Thus we obtain in the first order 



f (1) (x,ft) = exp 
whereas the second-order expansion reads 



1, {Am 

4 l0g Ur- 

9h 2 





uj 2 MQ 


-I 


8£l 2 Ah 




2 1 4 y 


4MQ 2 



UJ , , 



(3.5) 



■$ (2> (x,n) = exp 



1, fMtt\ 1 w 2 

T lOg — - — + -TTZ ~ 



UJ 



m 2 



\ hit J 

5 
2 



^2 r 



2r 



2057T 



1G 
3a; 2 
2fi 2 
2in 3 



4fi 2 

3?i 



Mfi /3 



4Mft 2 



3oj 2 uj a 
h 4ft 2 _ 8ft 4 

40 I 2 



2ft 2 



32Af 4 ft 6 2M 3 ft 5 ' 



llh 2 

8M 2 TF ; 



6Mft 3 



(3.6) 



Both in first and in second order, the ground-state wave function depends on the artificially introduced frequency 
parameter ft. According to the principle of minimal sensitivity ^] we minimize its influence on <]/( n ) (x. ft) by searching 
for local extrema of ^"^(a;, ft) with respect to ft. As we have written the wave function in the form ^f^(x, ft) — 
explW-™'^, ft)], it is sufficient to take into account just the inner derivative of \?W(a:, O) , i.e. we obtain the condition 
dW^(x,Q)/dn = 0. 

It turns out in the first order n = 1 that this equation has two solutions for x < 0.684 and for x > 0.780, however 
in the interval 0.684 < x < 0.780, ^^(x,fl) does not have any extremum Q. In accordance with the principle of 
minimal sensitivity we look for turning points on that interval instead, i.e. we solve d 2 W^{x,n)/dn 2 = 0. Fig. |l| 
shows how the curve for the turning points links the extremal branches. Now we have to choose which one of the 
branches of ft 1 ^) we take into account. Inserting the lower branch for x > 0.780 into the wave function (|3.5| ) leads 
to unphysical results as the ground-state wave function explodes dramatically. Thus we choose the upper branch for 
x > 0.780. For x < 0.684 the wave function becomes rather independent of the choice of ft. As we are looking for a 
function il^\x) which is as smooth as possible, we choose the lower branch for x < 0.684. Fig. |l| shows all branches 
of ft and highlights our final choice by a solid line. 
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FIG. 1. First-order results for the variational parameter £1 at the intermediate coupling g — 1/2. The extremal branches 
for x < 0.684 and for x > 0.780 (solid lines and dashed lines) are obtained from the equation dW^'(x, Q)/dQ = 0. For 

0. 684 < x < 0.780 there are no real positive solutions of this equation. Thus we look in this interval also for turning points, 

1. e. we determine real positive solutions of the equation d 2 W ( - 1 \x, Q)/dQ? = 0. The curve for the turning points on the entire 
interval lies between the two other branches (dot-dashed line) and fills the gap. Thus we can take those branches into account 
which provide us with the most continuous function Q^^x), i.e. the solid line. 



n (2, (x) 2 




FIG. 2. The two positive branches of f2^ 2 ' (x) on the interval [0,4] for intermediate coupling g = 1/2 obtained by solving 
the turning point equation d 2 W^ (x,0,)/dQ 2 = 0. In order to achieve the smoothest function we choose the lower branch for 
x < 0.8 and the upper branch for x > 0.8. This choice is justified by the results of our first-order calculation (see Fig. [j]). 

For the second order n = 2 there are no real positive solutions of the equation dW^ 2 '(x, fl)/d£l = on the interval 
x = [0,4]. Once more we have to look for turning points instead and solve d 2 W^ (x,Q)/dfl 2 — 0. This equation has 
two positive solutions on this interval, so we get two branches for the solution fl^(x) (see Fig. |J). Again we have to 
choose one of these two branches. Relying on a similar argument as for the first-order we choose the upper branch 
for x > 0.8 and the lower one for x < 0.8. The perturbation series converges so quickly that the curves for the first 
and second order as well as the exact ground-state wave function are not distinguishable on the plots. To see the 
difference we determine the mean square deviation from the exact numerical solution 



dx 



<$r(™)( 



* ex (a;) 



(3.7) 



where the index n denotes the order. The integration is performed numerically. The factor 2 is introduced for 
symmetry reasons, since we restrict our calculations to the positive x-axis. It turns out that the mean square 
deviation = 6.8 x 10~ 7 is smaller than = 1.1 x 10~ 5 by a factor of 0.063, which indicates that variational 
perturbation theory converges very quickly also for the ground-state wave function. The same applies to both weak 
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FIG. 3. The normalized second-order ground-state wave function for weak coupling (dashed, g — 0.1), for intermediate 
coupling (solid, g = 1/2), and for strong coupling (dotted, g — 50). 

and strong coupling as is illustrated in Fig. |^ which shows the second-order ground-state wave function ^^(x) for 
g = 0.1, g = 1/2, and g — 50, respectively. 

Note that variational perturbation theo ry do es not preserve the normalization of the wave function. Although the 



perturbative ground-state wave function (2.32) is still normalized in the usual sense, this normalization is spoilt by 
cxtrcmizing ^(x,fl) with respect to the frequency parameter f2. Thus we have to normalize the variational ground- 
state wave function at the end. 
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